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ABSTRACT 

We have used a new plasma simulation model to study 
the resonant interactions between streaming cosmic-ray 
ions and a self-consistent spectrum of Alfven waves, such 
as might exist in the interstellar medium upstream of a 
supernova remnant shock wave. The computational 
model is a hybrid one, in which the background 
interstellar medium is an MHD fluid and the cosmic- 
rays are discrete kinetic particles. The particle sources 
for the electromagnetic fields are obtained by averaging 
over the fast cyclotron motions. When the perturbed 
magnetic field is larger than about 10% of the 
background field, the macro- and micro-physics are no 
longer correctly predicted by quasi-linear theory. The 
particles are trapped by the waves and show sharp jumps 
in their pitch-angles relative to the background magnetic 
field, and the effective ninety-degree scattering time for 
diffusion parallel to the background magnetic field is 
reduced to between 5 and 30 cyclotron periods. We use 
our simulation results to suggest that Type I supernova 
remnants may be the principal sites of cosmic ray 
acceleration. 

Keywords: cosmic-ray acceleration, first-order Fermi 
acceleration, shock acceleration, upstream waves, 
quasilinear diffusion, supernova remnants 


1. INTRODUCTION 


In the late 1970's, several authors almost simultaneously 
suggested that cosmic-rays can be accelerated by a first- 
order Fermi mechanism, brought about as they scatter 
many times across astrophysical shock waves (Ref. 1). In 
the context of explaining the energy density and 
spectrum of galactic cosmic rays, the shock waves driven 
by expanding supernova remnants give the most likely 
setting for such acceleration. This hypothesis is 
attractive for at least two reasons: 1) it predicts a power- 
law particle spectrum similar to that deduced as the 
source spectrum for galactic cosmic rays, and 2) the 
overall energetics of galactic supemovae is sufficient to 
account for the observed cosmic-ray energy density, 


provided that the acceleration process is at least a few 
percent efficient. 

However the fact that the cosmic-ray source spectrum 
seems to extend smoothly up to 10 14 - 10 15 eV poses a 
potential difficulty for these theories. The reason is as 
follows. In diffusive cosmic-ray acceleration, a particle 
gains a small amount of energy each time it is scattered 
across the shock front. The highest-energy particles 
must cross and re-cross the shock front many times in 
the process of being accelerated. Thus considerable 
amounts of time are required to accelerate the highest- 
energy cosmic-rays. 

The effective upper limit on the energy to which a 
particle can be accelerated in these first-order Fermi- 
accele ration theories is given by equating the 
acceleration time to an estimate of the lifetime of the 
supernova remnant (Ref. 2). Using a minimum value of 
the spatial diffusion coefficient 

D mm = r L v / 3 - (1) 

where r L is the ion Larmor radius based on the 
background field, Lagage and Cesarsky 2 obtained an 
upper limit on the cosmic-ray energy: 

E max < a few X 10 13 (B 0 / 10" 6 gauss) Z/A eV. (2) 

Here Z/A is the nuclear charge-to-mass ratio of the 
cosmic-ray ion and B Q is the value of the interstellar 
magnetic field. This result is puzzling, since the cosmic- 
ray source spectrum is observed to be smooth in slope 
and continuous in value to energies in excess of this 
limit, suggesting that there is not a transition in 
acceleration mechanisms at the energy given by equation 
(2). One would like to establish a more realistic 
maximum cosmic-ray energy by using a self-consistently 
derived diffusion coefficient to replace that in Eqn. (1). 
The computational simulations we describe here can 
begin to address this issue by explicitly evaluating the 
diffusion coefficient of fast particles in a self-consistent 
spectrum of Alfven waves, and studying its scaling with 
cosmic-ray flux. 
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The efficiency of diffusive cosmic-ray acceleration is also 
an important issue. First-order Fermi acceleration 
predicts that the efficiency depends on the spatial 
diffusion coefficient, which itself depends on the 
scattering rate. Frequently, quasi-linear assumptions 
have been used to describe scattering in the self- 
consistent spectrum of Alfven waves upstream of the 
shock front. Yet several different arguments (Ref. 3) 
suggest that the Alfven-wave perturbations near the 
shock are large (S B/B Q > 1), thus violating the initial 
quasi-linear assumptions. Recent observations of fast 
particles and Alfven waves at interplanetary shocks (Ref. 
4) show several types of non-quasi-linear behavior. The 
computer simulation study described here allows a more 
general examination of these effects, since it does not 
require quasi-linear conditions to be satisfied. 


2. THE COMPUTATIONAL MODEL 

The intention is to model a region upstream of a 
supernova remnant shock-wave, in which there exists an 
anisotropic distribution of cosmic-rays as a consequence 
of the previous action of the shock-acceleration process. 
The shock-wave itself is not included in the simulation 
volume; its presence is reflected via the anisotropy of the 
initial cosmic-ray distribution. We study the time 
evolution of this initially anisotropic cosmic-ray 
distribution, together with the self-consistent Alfven 
waves which it produces. 

Our simulation model is an orbit-averaged Darwin 
quasi-neutral hybrid code (Ref. 5). In less technical 
terms, this means the following. The electrons and ions 
in the interstellar medium are treated as a single ideal 
MHD (magnetohydrodynamic) fluid, while the cosmic- 
ray ions are treated as discrete kinetic particles. The 
displacement current is assumed small relative to plasma 
currents, and the cosmic-ray number density is assumed 
small relative to the density of the interstellar medium. 

The simulation geometry is one-dimensional, with a 
periodic grid having either 128 or 256 mesh points. 
Magnetic and electric fields are defined on the grid 
points. The direction of spatial variation (the x-axis) is 
parallel to the background magnetic field, so that all 
waves propagate parallel or antiparallel to Bq. All 
perturbed magnetic fields are perpendicular to the 
background field. 

Cosmic-rays are represented by 32000 discrete particles, 
whose currents contribute as sources for the self- 
consistent magnetic and electric fields at the grid points. 

Particle positions are defined along the whole x-axis, not 
just at grid points. Particle trajectories in position and 
velocity are integrated via micro-timesteps At, with 
resolution sufficient to follow the cyclotron motion in 
detail. Currents due to the particles are averaged over 
much larger time intervals aT, where typically N - 
AT/At = 88. This orbit-averaging method allows the 
field-solve to be done much less frequently than the 
particle-push. In addition, it results in a discrete-particle 
noise level that is lower by a factor of N than that which 
would otherwise obtain for the same number of cosmic- 
ray particles in the absence of orbit-averaging (Ref. 6). 
As a consequence, it is practical to simulate instabilities 
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Figure 1: Initial cosmic-ray distribution function for 
Models A and B. Contours show equal values of p^Cp) ; 
vertical axis is p/mc = 7 v/c ; horizontal axis is the cosine 
of the particle's pitch-angle relative to Bq. 

which grow slowly relative to the cyclotron period of the 
cosmic-rays. 

The initial distribution function of the cosmic-rays is a 
power-law in the magnitude of the momentum, with 
spectral index 4: 

f(p) d 3 p p -4 d 3 p , 0.1 < 7 v/c < 2.0 . (3) 

This power-law momentum distribution is then taken to 
be drifting relative to the laboratory frame with a (non- 
relativistic) drift velocity v D = 0.1 c = 10 v A , where v A 
is the Alfven speed: 

v A 2 = B 0 2 / 4 * p . ( 4 ) 

Figure 1 illustrates the initial shape of this drifting 
power-law cosmic-ray distribution function in our 
simulations. 

We shall discuss in some detail two contrasting 

simulations, Model A and Model which differ only in 
the initial density of cosmic-rays relative to the 

background interstellar medium: 

Model A: n cr /n Um = 4 x 1(T 3 ; n Cf v D / n ism v A = 0.04 ; 

Model B: njn^ = 4 x 10’ 2 ; n cr v D / n ism v A = 0.40 (5) 


3. RESULTS OF THE SIMULATIONS 


The overall evolution in both Model A and Model B is 
as follows. The initial anisotropy of the cosmic-rays 
relative to B Q drives a resonant Alfven-wave instability in 
the interstellar medium. A broad spectrum of Alfven 
waves grows up in time, since modes with different 
wavenumbers k can resonate with different particles in 
the power-law distribution function. When the Alfven 
waves reach large amplitude, mode-coupling between 
them also drives sound-wave fluctuations in the 
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Figure 2: Time histories of energy density in the most unstable Alfven-wave mode, 
a) Model A, n^/n^ = 4 x 10' 3 ; b) Model B, n cr /n jsm = 4 x 10' 2 . Vertical axis is 
kB k 2 /B 0 2 , horizontal axis is cyclotron frequency times time. 


background interstellar medium. As a consequence of 
the Alfven-wave fluctuations, the mean cosmic-ray drift 
velocity decreases from its initial level as individual 
cosmic-rays are scattered in pitch-angle. 

Frames a and b of Figure 2 show the time evolution of 
energy in the largest unstable Alfven-wave mode for 
Models A and B respectively. An initial period of 
exponential growth of the resonant anisotropy-driven 
instability is followed by saturation of the Alfven-wave 
amplitudes. For Model A (the low-flux case), root- 
mean-square magnetic field perturbations in the 
saturated phase are (6 B/B 0 ) rms * 3 x 10' 2 ; individual 
magnetic-field peaks with amplitudes 0> B / B o)peak “ 
0.15 are present in snapshots of 6 B/B Q as a function of 
x. For Model B (the high-flux case), rms magnetic field 
perturbations in the saturated phase are (6 B/B 0 ) rms - 
0.2 - 0.3 , with individual peak field amplitudes 
(sB/Bq)^^- 0.9 - 1.1 . Thus one interesting lesson is 
that the description of a spectrum of Alfven-waves in 
terms of rms perturbation levels can lead one to forget 
that individual magnetic-field maxima are present at 
values many times larger than the rms level. We shall 
see below that these larger-than-average peaks can have 
a strong influence on the non-quasilinear behavior of the 
particle scattering. 

For Model A, the low-flux simulation, the Alfven-waves 
propagate forward (in the direction of the cosmic-ray 
anisotropy). However for Model B, where the Alfven- 
wave amplitudes are larger, we observe some backward 
propagating Alfven-waves. We also observe sound-wave 
perturbations in the background interstellar medium, 
driven by wave-wave coupling between different Alfven 
wave modes. We have suggested elsewhere (Ref. 7) that 


sound-waves driven by Alfven-wave turbulence upstream 
of astrophysical shocks may contribute to the electron 
density fluctuations inferred from interstellar scintilla- 
tion measurements. Figure 3 illustrates the Alfven-wave 
amplitudes (left panel) and the sound-wave amplitudes 
(right panel) as functions of x and t, for Model B. 
Forward propagation corresponds to a line moving 
upwards and to the right in these pictures. 

As a result of scattering of individual cosmic-rays by the 
Alfven waves, the mean drift velocity v D of the cosmic- 


x 

Figure 3: Amplitudes of Alfven waves (left panel) and 
sound waves (right panel) in Model B, as functions of x 
and t. 
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Figure 4: Time evolution of cosmic-ray drift velocity, 
v D , in Model A. Horizontal axis is cyclotron frequency 
times the time. 


rays decreases with time. Figure 4 shows the decline of 
Vp/c with time for Model A, the low-flux simulation. 
Tne horizontal axis ft c t is the cyclotron frequency times 
the time. 


the type of behavior that is predicted by quasi-linear 
theory. Note, however, that the particle diffuses through 
/i = 0 without any hesitation, whereas lowest-order 
quasilinear theory predicts that diffusion through ^ = 0 
should not be possible. 

The test-particle in Model B, the high-flux case, shows 
markedly different behavior. Its trajectory is dominated 
by large jumps in pitch-angle, S/Vm ~ 1 , within 5 to 10 
cyclotron periods, and by repeated passes through n = 
0 . These large jumps in pitch-angle are decidedly non- 
quasilinear, since they involve order-unity changes in 
pitch-angle in a few cyclotron periods. To investigate 
their physical origin, we followed in detail all of the 
forces experienced by a test-particle during one jump 
from n « + 1 to /i ~ -1. The largest systematic force is 
the x-component of the Lorentz force, ( 6v± x 5B± ) x , 
representing the nonlinear force due to the perturbed 
particle velocity and the perturbed magnetic field. 


We have used test-particles to study the scattering 
process of cosmic-rays in these self-consistent Alfven- 
wave spectra. In frames a and b of Figure 5, we plot the 
trajectories of two typical test-particles in Models A and 
B respectively. The vertical axis ^ in these plots is the 
cosine of the particle's pitch-angle relative to the actual 
magnetic field; the horizontal axis nt is the cyclotron 
frequency times the time. We have constructed 
analogous trajectory plots using the pitch-angle relative 
to the background magnetic field B Q , which is more 
relevant to calculation of the spatial diffusion 
coefficient; we find similar results to those shown here. 

The test-particle in Model A, the low-flux case, shows a 
trajectory that is qualitatively diffusive in nature, made 
up of a sequence of small changes in pitch-angle. This is 


Figure 6 shows the time evolution of this ( 8 vx x & Bx ) x 
force on the test-particle during its rapid jump from p ~ 
+ 1 to /i ~ -1 . Note that the timescale on the horizontal 
axis is greatly expanded relative to that in Figure 5. 
During the rapid jump, the particle experiences a series 
of brief intervals when the ( S\± x 6 B± ) ? force is large, 
both in the positive and negative directions. However 
during this particular time interval there are four large 
negative spikes in («vix fiBi ) x which outweigh the 
smaller positive spikes, and which together are 
responsible for the large net force on the test particle in 
the negative x direction. One interpretation of these 
large negative spikes is that they represent intervals 
during which the test-particle experiences magnetic 
trapping in large-amplitude Alfven waves. This is 
supported by the analysis in Ref. 6 of the phase angle 
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Figure 5: Trajectories of test particles showing evolution of pitch-angle p = cos 0 , as 
function of the cyclotron frequency times the time, a) Test particle in Model A (low 
cosmic-ray flux); b) The same test particle in Model B (high cosmic-ray flux). 
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Figure 6: The x-component of the v x B force on a test 
particle in Model B, during the brief time interval when 
the particle is undergoing a sharp negative jump in pitch- 
angle. Vertical axis is x-component of the Lorentz force, 
normalized to B 0 ; horizontal axis is cyclotron frequency 
times the time, using a greatly expanded scale. 


between v± and B± , which pauses near zero degrees 
during the main negative spikes shown in Figure 6. 

We conclude from these results that the higher-flux 
simulation, Model B, shows substantial non-quasilinear 
behavior, particularly in its particle scattering 
characteristics. To assess the effect on the spatial 
diffusion coefficient of cosmic-rays, we have computed 
the effective ninety-degree scattering time in our 
simulations using an ensemble of test-particles. We 
compute the ninety-degree scattering time as a function 
of the magnitude of the particle momentum, for each of 
the simulation conditions we have modeled. The 
calculation is based on constructing the auto-correlation 
function of the pitch-angle relative to B Q , followed along 
a test-particle's trajectory: 

T 

max 

c„( o = (1 /tj ; o*(t)- <M>)o*(t+o- </*>)dt (6) 

0 

The autocorrelation time r ac is calculated from the 
half-width of the central peak of C^( r) . The ninety- 
degreescatteringtime is then given by 

-90 “ < -ac /C M < r= °) ( 7 > 

where the brackets indicate an ensemble average over 
all of the test-particles. Further details of this 
calculation, as well as graphs showing the 
autocorrelation functions themselves, may be found in 
Ref. 6. 


4. IMPLICATIONS FOR COSMIC-RAY 
ACCELERATION 


4.1 Scaling Laws 

Computational realities have constrained us to simulate 
only a relatively small portion of the parameter space of 
interest to cosmic-ray acceleration at supernova remnant 
shock-waves. Thus, for example, our results do not apply 
directly to the scattering of the highest-energy cosmic- 
rays, since the simulations lack sufficient dynamic range 
in particle energy. However, we can use our simulation 


results to construct scaling laws and upper limits for the 
important physical quantities involved in diffusive 
acceleration. These scaling-laws give a suggestion of the 
physical behavior that will apply for parameter ranges 
which are presently out of reach of direct computational 
simulation. 

In particular, the values of the ninety degree scattering 
time constructed from our simulations can be used to 
investigate the scaling of the cosmic-ray diffusion 
coefficient. The spatial diffusion coefficient is given by 

0 -^/ 3 ; ( 8 ) 

we can use the values of from our simulations to 
investigate D . First we consider the variation of r ^ 
with the amplitude of the magnetic-field perturbations. 
Figure 7 shows the dependence of the cyclotron 
frequency times the ninety degree scattering time, fir™ , 
on the root-mean-square strength of the magnetic field 
perturbations, denoted here by sB/B 0 . The three 
curves shown in Figure 7 correspond to three different 
values of the magnitude of the particle momentum, 
7 v/c . An approximate fit to these curves is given by the 
following relations: 

fir 90 - K i( * B /B 0 )' 3 ±0 ' 6 , 0.4 < s B/B 0 < 0.7 ; 

Or go . Kj ( S B/B 0 )-°- 4 , 0.06 < S B/B 0 < 0.4 ; (9) 

where 1.1 x 10 1 102 , and 140. 

Recall that the quasi-linear model normally used for 
calculations of cosmic-ray scattering is that of Jokipii 
(Ref. 8), which predicts that 

fir 90“ (4A) (B 0 2 /kB k 2 ), (10) 

where B k is the Fourier amplitude of the magnetic field 



logio(<5B/B 0 ) 

Figure 7: Scaling relation between ninety degree scat- 
tering time and root-mean-square amplitude of mag- 
netic field perturbations $B/B 0 , from a collection of 
simulations. Lines represent ninety-degree scattering for 
different particle momenta: p/mc = 0.1 (filled tri- 
angles), 02 (filled squares), and 0.4 (empty triangles). 
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fluctuations at wavenumber k. If we roughly equate 
IcB^/Bq 2 in Eqa (10) with («B/B 0 ) 2 , then Jokipii's 
quasilinear theory would predict n « ($B/B 0 )' 2 in 
Eqn. (9), which should apply at the lower values of 
fiB/B 0 . Instead, the observed dependence a 

(sB/Bq)- 0 - 4 for small values of 6 B/B 0 is considerably 
weaker than this prediction. We conclude that for the 
perturbation amplitudes we have studied ( 6 B/B Q > 
0.06 ) , the ninety-degree scattering times seen in our 
simulations do not follow the simple quasilinear trends 
generally assumed in theories of cosmic-ray diffusion. 

4.2 Maximum Energy of C osmic-Ravs 

As a substitute for quasi-linear theory, we can use the 
empirical relations given in Eqn. (9) and Figure 7 to 
calculate E , the maximum cosmic-ray energy that 
can be produced during the active lifetime of a 
supernova remnant. We shall follow the general 
approach of Lagage and Cesarsky (Ref. 2), using our 
empirically derived scaling for the ninety degree 
scattering time. We calculate the cosmic-ray momentum 
as a function of time, using the Sedov solution for the 
evolution of the supernova remnant shock-wave, and 

find the asymptotic maximum cosmic-ray energy E max 
together with the time needed to achieve this energy. 

Since we are not able to simulate directly the scattering 
of particles with energies typical of the highest-energy 
cosmic rays ( 10 14 - 10 15 eV ), we shall have to make 
some assumptions concerning the self-consistent 
behavior of these particles and of the Alfven waves 
which scatter them. According to resonant quasi-linear 
theory, particles at a given momentum p generate the 
waves by which they are themselves scattered, at the 
resonant wavenumber k(p) = 1 / r L (p) « 7" 1 • The 
highest-energy cosmic-rays are rare compared with those 
at lower energies, because of the power-law distribution 
function (Eqn. 3). As a result, one might expect the 
wave amplitudes at small k (long wavelength) to be 
smaller than those at higher k. To derive an up per lim it 
for E max , we neglect this wavelength distinction and 
assume that the scaling results of our simulations also 
apply to particles having higher energies than those 
present in the simulations. 

Assuming a spatially uniform diffusion coefficient, 
Lagage and Cesarsky found the rate of gain of 
momentum for a particle encountering a shock in a 
diffusive manner to be 

p' 1 dp/dt = ( Uj - u 2 )( Uj u 2 )/[ 3D ( u x + u 2 ) ] (11) 

where u x and 02 are the fluid velocities up and 
downstream of the shock, measured in the shock frame. 
We assume a strong shock, so that u 2 /uj = 1/4 . If we 
substitute our results from Eqn. (9) into Eqn. (8) for D, 
and then use this diffusion coefficient in Eqn. (11), we 
find for relativistic particles that 

dp/dt = ( 3/20 ) ( eB 0 /nr^ ) ( Ul / c ) 2 

= (3eB 0 /20K u )(<B/B 0 )“(u 1 /c) 2 (12) 


where Kj and Kj are defined below Eqn. (9) and 
a = 3± 0.6 , 0.4< «B/B 0 < 0.7 , 
a =0.4, 0.06 < «B/B q < 0.4 . (13) 

Here B 0 is the magnitude of the background magnetic 
field upstream of the shock. 

We evaluate u x in Eqn. (12) using the Sedov solution for 
the expansion history of the supernova remnant. This 
means we assume in this discussion that the cosmic-ray 
pressure is less than the ram pressure at the shock front, 
so that the shock speed is not strongly modified by the 
presence of the cosmic rays. If we normalize to 
parameters typical of Type II galactic supernovae and of 
a hot Phase III interstellar medium, then the upstream 
velocity in the shock frame is 

Ul - (10 3 5 km s' 1 ) x (Vgj/10 3-5 km s' 1 ) (t/t^ 3 / 5 (14) 

where v • is the initial ejecta expansion velocity and t t is 
the timetor the ejecta to sweep up their own mass from 
the interstellar medium: 

t t « (10 3 5 years)(E SN /10 51 ergs) 1 / 3 (n Um /10' 2 cm’ 3 )' 1 / 3 

x (Vgj/10 3,5 km s" 1 )' 5 / 3 . (15) 

Substituting expression (14) for u, into Eqn. (12), we can 
integrate to find the time dependence of the cosmic-ray 
momentum: 

p(t) = R. [ 1 - (0.91 t x / t) 1 / 5 ] (16) 

where 

p M = (7 x 10 12 eV/c)x (B 0 /10' 6 gauss) (11Z/AK 12 ) 
x (s B/B Q ) “ (E^/IO 51 ergs) 1 / 3 (n ism /10- 2 cm- 3 )' 1 / 3 
x (Vej/10 3 * 5 km s' 1 ) 1 / 3 (17) 

According to Eqn. (16), a cosmic-ray reaches a 
momentum equal to 50% of at a time t = 29 t j * Kr 
yrs, for the nominal parameters used to normalize the 
above expressions. However the asymptotic momentum 
attainable, p,,,, will be considerably less than the 
observed value of 10 14 - lO 1 ^ eV/c unless 6 B / B« > 1 
and some of the physical parameters in Eqn. (17) differ 
substantially from the nominal values to which they have 
been normalized here. 

The parameter $B/Bp is of course a function of the 
actual physical conditions near the supernova remnant 
shock, and must be determined self-consistently as a 
function of the externally given environment. We can 
use our simulation results to relate 6B/B 0 to the up- 
stream flux of cosmic-rays, J^, defined dimensionlessly 
as 

J cr" ( n cr v cr)/( n ism v A)' ( 18 > 

Figure 8 shows the dependence of S B/B„ on J CT in our 
simulations. The result is a strikingly simple power-law, 
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Figure 8: Scaling relation betwen root-mean-square am- 
plitude of magnetic field perturbations $B/B 0 , and 
normalized cosmic-ray flux , from our simulations. 


S B/B Q = 0.57 J CT 0 ' 7 • (19) 

If we substitute this relation into Eqn. (8) for the cosmic- 
ray diffusion coefficient, we find that 

D = (r L v/3)(59J cr - 21 ) (20) 

for 1.3 > J cr > 0.6 . If one extrapolates this relation to 
somewhat larger fluxes, we obtain the result that 
D < r L v /3 for J cr > 7 . Previous authors (e.g. Refs. 2 
and 3 and references therein) have suggested that D = 
r L v /3 is a minimum value for the cosmic-ray diffusion 
coefficient. However the mean free path can be smaller 
than the Larmor radius in the background field, r, , 
because the fluctuating fields can be larger in magnitude 
than the background field when J cr > 3, and because 
transitory strong trapping in the turbulent magnetic 
fields can reduce the mean free path while still allowing 
the particles to diffuse. 

We can substitute Eqn. (19) into expression (17) for p^ , 
and derive the asymptotic maximum cosmic-ray energy 
as a function of J cr . For simplicity we take Z/A = 1, 
since the mechanism must work for protons. Then in the 
range 0.6 < J cr < 1.3 , we have 

p„ = (1.4 x 10 12 eV/c)x J CT 21±a4 (B 0 / 10- 6 gauss) 

x (E sn / 10 51 ergs) 1 / 3 (n jsm /10' 2 cm" 3 )- 1 / 3 

x (Vgj/10 3 - 5 km s' 1 ) 1 / 3 ; (21) 

while for 0.04 < J cr < 0.6 , we have 

p„ = (0.45 x 10 12 eV/c) x J^ 0 - 28 (B Q / 10- 6 gauss) 

x (Esn/ 10 51 ergs) 1 / 3 (n^/10' 2 cm’ 3 )’ 1 / 3 

x (Vgj/10 3 - 5 km s' 1 ) 1 / 3 . (22) 

The numerical value of the asymptotic cosmic-ray 
momentum p^ predicted by Eqns. (21) and (22) depends 
on the cosmic-ray flux J CT , and weakly on the physical 


characteristics of the supernova remnant and its 
surroundings. First consider the variation of p,, with , 
and assume that the other physical parameters take on 
the normalization values indicated in Eqns. (21) and 
(22). Then over the range we have simulated, 0.04 ^ J 
< 1.3 , the factors containing in Eqns. (21) and (22) 
lead to a relatively narrow variation of p^ : 0.2 x 10 12 < 
p^ ^ 2.4 x 10 12 eV/c . The observed maximum cosmic- 
ray momentum is considerably larger than this range. In 
order to produce realistically large cosmic-ray momenta 
by varying J CT alone, one would have to extrapolate our 
results to the regime J CT > > 1, where 5 B/Bq > 1. 

Because of the weak dependence of p^ in (21) and (22) 
on the parameters of the supernova explosion and the 
interstellar medium, J CT > 8 is the only way consistent 
with the normal energetics and velocities of supernovae 
to raise p ro to its observed value ~ 10 15 eV. Under such 
circumstances the cosmic-ray pressure can alter both the 
shock structure and the value of u t used in Eqn. (14). 
For the p“ 4 momentum spectrum expected in the 
simplest form of shock acceleration, with the minimum 
momentum p . comparable to the thermal momentum 
of the shockea interstellar medium, p cr < pv 2 implies 

J cr < 0.9 (100 keV/E min ) 1 / 2 (v^/10 3 - 5 km s' 1 ) 3 

x (n^/10' 2 cm' 3 ) 1 / 2 (B 0 /10' 6 gauss)' 1 / 2 

x { 33/ [ 1 + 21n(p max /rnc) ] } (1,/t) 6 / 5 (23) 

Here > E min = Pmin 2 /2m is the kinetic energy 
corresponding to the lower cutoff of the momentum 
distribution. The fact that x-ray emission from galactic 
supernova remnants is in rough accord with the 
predictions of supernova remnant models having p CT < 
pv 2 upstream (Ref. 3) therefore suggests that J > > 1 is 
unlikely in Type II supernovae. However, the higher 
velocity Type I supernova remnants (v e j = 10 4 km s' 1 ) 
have a much larger kinetic energy flux upstream from 
the shock for the same age relative to t,, and thus may 
be able to have high cosmic ray flux while still satisfying 
p cr < pv 2 ; J CT as large as 25 is allowed by Eqn. (23), 
corresponding to 6 B/B^ up to 5.5 in the cosmic ray 
foreshock. From Eqn. (21), we find p^ ~ 10 15 eV/c if J 
~ 19, corresponding to very strong turbulence 6 B/B Q « £ 
a value marginally within the dynamical constraints. 
Thus, on dynamical grounds, Type I supernova remnants 
may be the most favorable site for shock acceleration. 
One should note that observations of the rate of 
supernova occurence, lifetime of cosmic rays in the 
galaxy inferred from the abundances of heavy nuclei in 
the cosmic rays, and the number density of cosmic rays, 
can be used to suggest that < 1 even in Type I 
remnants during the acceleration phase (Zachary, Max 
and Arons, in preparation). This conclusion militates 
against the formation of large $B/B 0 even in Type I 
supernova remnants. However, recent work (e.g.. Ref. 
9) has suggested that the protons may have undergone a 
different history in the galaxy from that experienced by 
the heavy nuclei. If the proton lifetime is a factor of ten 
shorter than that of the heavy nuclei, our results suggest 
that shock acceleration in Type I remnants remains a 
viable candidate for the origin of cosmic rays. 

In summary, extrapolation of our simulation results to 
high-energy cosmic rays implies an upper limit of about 




52 


C E MAX &AL 


10 12 eV for the maximum energy attainable in shock- 
acceleration at Type II supernova remnant shocks. This 
is several orders of magnitude less than the observed 
maximum energy of the power-law source distribution of 
cosmic rays. On the other hand, the expansion of Type I 
supernovae into the low density phase of the interstellar 
medium may be able to accelerate cosmic rays to 
energies comparable with the observed upper cutoff of 
the smooth power-law distribution, if the very large field 
amplitudes implied by our simulations are applicable to 
the long wavelength magnetic fluctuations with which 
the highest energy particles interact. 
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